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Abstract 

We present both theory and an algorithm for solving time-harmonic wave prob- 
lems in a general setting. The time-harmonic solutions will be achieved by com- 
puting time-periodic solutions of the original wave equations. Thus, an exact con- 
trollability technique is proposed to solve the time-dependent wave equations. We 
discuss a first order Maxwell type system, which will be formulated in the frame- 
work of alternating dilfcrcntial forms. This enables us to investigate dilfcrcnt kinds 



o 

(N 

< 

^ of classical wave problems in one fell swoop, such as acoustic, electro-magnetic or 

^ tt^ elastic wave problems. After a sufficient theory is established, we formulate our ex- 

act controllability problem and suggest a least-squares optimization procedure for 
^ its solution, which itself is solved in a natural way by a conjugate gradient algo- 

^ rithm operating in a purely L^-type Hilbert space. Therefore, it might be one of the 

0^ biggest advances of this approach that the proposed conjugate gradient algorithm 

does not need any preconditioning. 

Key Words wave equation. Maxwell's equations, differential forms, differential 
geometry, time-periodic waves, time-harmonic waves, controllability, least-squares 
formulation, conjugate gradient method, discrete exterior calculus, discrete differ- 

ential forms 
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1 Introduction 

Time-harmonic wave propagation is an important phenomenon which has many ob- 
vious applications in acoustics, electro-magnetics and elasticity, among others. Tradi- 
tionally, the numerical solution approaches have been based on finite differences, finite 
elements or boundary element techniques. As our goal is to consider heterogeneous media 
as well, we pay attention to methods based on partial differential equations. Hence, some 
kind of tessellation of the spatial domain is necessary. 

To obtain accurate results for wave propagation, the discretization mesh needs to be 
adjusted to the wavelength. If the time-harmonic case is directly addressed, one is faced 
with the solution of a large-scale indefinite linear system which is a difficult task. 

Instead of solving directly the time-harmonic problem for a given frequency u G IR+, 
it is possible to compute the solution by control techniques. Then the solution is found 
by searching for an appropriate initial data for the wave equation which minimizes a 
quadratic functional that measures the difference between the initial state and the final 
state after one time period T = 2'7i/u. A natural quadratic error functional is the squared 
energy norm of the system, allowing to minimize the cost by the conjugate gradient 
method (CGM) operating in Hilbert spaces. This approach has been successfully applied 
to acoustics, electro-magnetics and elasticity [IUl[IIllI21[I71[IHl[IHl[2ni[2Il[221[231[23l2Sl 
[26] . In practice, the method seems to have a good asymptotic computational cost. Even 
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though no theory exists, the computational cost of the method seems to be of order 0(n), 
where n is the number of spatial degrees of freedom. The drawback of using the traditional 
(second order in time) formulation of the wave equations is that the energy norm is then 
of H^-type, and as such, the minimization problem is badly conditioned. This is handled 
by applying preconditioning to the conjugate gradient minimization. Unfortunately, this 
means that a discrete elliptic problem (linear system) still needs to be solved at every 
conjugate gradient iteration step. In recent papers, the linear system has been solved by 
an algebraic multi-grid method which still maintains the good asymptotic performance of 
the solution technique, but makes it quite more difficult to implement the solver to utilize 
the computing power of modern parallel computers and multi-core processors. 

Hence, an alternative approach has recently been proposed in the short paper of 
Glowinski et al. [22]. The idea is to formulate the control method for an equivalent 
first order system which has an L^-type energy norm, and hence, a well-conditioned min- 
imization problem results. This eliminates the need for preconditioning the conjugate 
gradient minimization and, thus, greatly simplifies the parallel implementation of the 
method. This approach also has drawbacks as the spatial discretization needs to be 
based, for example, on mixed finite elements like Raviart-Thomas elements, which are 
more difficult to implement than standard finite elements. Initial numerical experiments 
(still unpublished) support the hypothesis that the cost of the new approach is also of 
order 0(n). 

In our project, we aim at generalizing the approach of [22] to generalized Maxwell 
equations formulated in terms of differential forms. The same formulation covers electro- 
magnetic, acoustic and elastic waves and it can be naturally discretized by so-called 
discrete differential forms (DDF) or discrete exterior calculus (DEC), which has recently 
been under very active research [221 [IS [IS] • Our goal is to develop theory and software for 
efficiently solving the generalized Maxwell equations using a control approach. We present 
a new solution theory for the generalized Cauchy problem (CP) at hand, such that we 
can be sure to have uniquely determined solutions evolving in time. Here the papers 
[121 US |35l [361 [33 [Ml [3l] as well as the monograph [31] are useful. Moreover, theoretical 
questions about the domain truncation procedure have to be considered. We are planing to 
use absorbing boundary conditions (ABC), generalized Dirichlet-to- Neumann operators 
(DtN), i.e. electric-to-magnetic operators (EtM), as well as perfectly matched layers 
(PML). All these techniques have to be developed for differential forms. The resulting 
software is targeted to mid-frequency variable coefficient wave propagation problems in 2D 
and 3D domains, where the dimension of the computational domain is 10-100 wavelengths. 
The software is targeted for modern parallel computers and multi-core processors. 

In this first paper we present and explain the basic ideas of our control approach for 
wave equations formulated as first order systems and using generalized Maxwell equations 
in terms of differential forms. First, in section |4] we investigate the Cauchy problem 
and establish a solution theory which meets our needs utilizing the spectral calculus for 
unbounded selfadjoint linear operators in Hilbert space. Then, in section [5] we introduce 
the least squares formulation and discuss the derivative of the least squares functional, 
which is the essential ingredient in our resulting algorithm, since we plan to use a conjugate 
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gradient method. In section[6| we discuss the conjugate gradient algorithm (CGA) in some 
detail. In section [7| we translate our results presented in terms of differential forms to the 
classical framework of vector analysis. We briefly demonstrate, which classical problems 
are covered by our generalized theory. Finally, in section |A] we present some preliminary 
numerical results and in section [8] we outline the ongoing work in this project. 

2 Notations and preliminaries 

We investigate wave scattering problems taking place in an exterior domain Q of the 
Euclidean space M^, which will be considered as a smooth A^-dimensional differentiable 
Riemannian manifold with a Lipschitz boundary T := dQ. 

o 

We define the space C°°''^{Q) of C°°-g-forms with compact support in Q. This space 
admits a natural scalar product 

{E,H)^ {E,H)n:= I EA*HeC, 

Jn 

where * denotes the Hodge star operator with respect to the Euclidean metric in M^, A 
the wedge product and the bar complex conjugation. Using this scalar product and its 

o 

induced norm we may define L^'''(fi) as the closure of C°°''^{Q). Then L^'''(f2) equipped 
with the scalar product ( ■ , ■ )L2,<i(n) '■={■, ■)n becomes a Hilbert space, the Hilbert space 
of square integrable g- forms on Q. 

As usual, we denote the exterior derivative by d and the co-derivative by 5. Thus we 
have on g-forms 

5= (-l)('?-i)^*d*. 

With respect to the latter scalar product the linear operators d and 6 are formally skew- 

o o 

adjoint to each other, i.e. for pairs of forms {E,H) G C°°'^(f2) x C°°'^"*"^(r2) we have by 
the weak version of Stokes' theorem 

0= [ d{EA*H)= [ dE A*H + {-ly [ EAd*H 
Jn Jn Jn 




= (dE,iy)L2,,+i(n) + {E,6H)i2,,^Qy 

This yields the possibility for weak versions of d and S (in the sense of L^(f2)-valued 
distributions) using smooth, compactly supported forms as test-forms. Hence, as usual 
we may define dE for a L^'''(f2)-form E and say E has weak exterior derivative, if there 

exists a l^^''+\n)-ioTm G, such that for all $ G C°°^i+\Q) 



Computation of Generalized Time-Periodic Waves 



5 



holds. Of course, we may define a weak co-derivative in the same way. Then we put 

D\Q) := {E e L2'«(Q) : dE e 

A^iQ) := [H e l^'^iQ) : 6H G 

Equipped with their natural graph-norms these are Hilbcrt spaces. Furthermore, we 
generalize the (electric) homogeneous boundary condition, which models a perfectly con- 
ducting obstacle and means that the tangential trace l*E of a differential form E vanishes, 
where t : T ^ Q denotes the natural embedding of the boundary manifold F regarded as 

o 

an {N — 1) -dimensional Riemannian submanifold of Q. For this purpose we define D^(Q) 

o 

to be the closure of C°°''(Q) in the norm of D^(Q). Indeed, by Stokes' theorem and a 
density argument one may easily check that for sufficiently smooth forms a vanishing tan- 

o o 

gential trace is generalized in D^(Q). D^(Q) is also a Hilbert space as a closed subspace of 

o 

D^(Q). An index at the lower left corners of the spaces D^(Q), D^(Q) or A^(Q) indicates 
vanishing exterior derivative or co-derivative, respectively. 

Another way to define these Hilbert spaces is to look at the densely defined linear 
operator 

do : £^'"(0) C L2'9(0) L2'9+i(0) 
E I — > dE 

and its adjoints, which will be marked by stars. Then do = d** is the weak exterior 

o o 

derivative on its domain of definition -D(do) = D^(Q). The kernel of do equals oD^(f2). Its 
adjoint operator d* = do equals by definition the negative weak co-derivative —6 on its 
domain of definition A^+^(f2), i.e. 

-dl^:S : A«+i(0) C L2'9+i(0) L2'«(0) 

H I — > 6H ' 

This is easy to see: Let H G D{d*) and d^H = F. Then by definition 

WE G L'(do) {dE,H)^2,,+i^a) = (£;,F)L2,<,(n), 

which is just the definition of the negative weak co-derivative. Therefore, H is an element 
of D{dl) = A'^+\n) and dlH = -5H holds. 
Since SS and dd vanish in the smooth case, 

dd = 0, 55^0 

still hold true in the weak sense and we also have the well known and important formula 

d5 + 5d = A, 

where the action of the Laplacian is to be understood componentwise with respect to 
Euchdean coordinates. Moreover, we get with closures taken in L^'''(Q) 

d D«-i(Q) c oD«(Q), SAi+^{n) c oA^(^^). 
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Let us formally define matrix-operators 
M : 



"0 


6 


, A:= 


'e 0" 


d 





/i 



Ma := iA-^M, 



where e respectively /i is a real, linear, symmetric, bounded and uniformly positive definite 
(with respect to the L^''^(f2)- respectively L^''^+^(i7)-scalar product) transformation on q- 
respectively {q + l)-forms, which is independent of time, and i denotes the imaginary 
unit, e and /i model material properties, i.e. in classical electro-magnetic theory e is the 
dielectricity and /i the permeability of the underlying medium. We note that e and /i 
are even allowed just to have L°°(i7)-entries in their matrix representations I'j, j given by 
chart bases 

uE = ^u^, jEjdhj, if E = ^Ejdhj. 

J 

Since do and 6 are skewadjoint to each other in this setting the unbounded linear 
operator 



Ma : D{Ma) C L^''''^+'(fi) 
iE,H) 



Ma{E,H) 



(2.1) 



with 



L^''^'''+'(r]) := L2''?'9+1(^]) := {^^^{n) X 
(as a set) equipped with the weighted scalar product 

where ( ■ , ■ )L2.'j.<7+i(n) denotes the canonical scalar product in the product space L^'''''^"'"^(ri), 
and domain of definition 

d(Ma) := h\n) X /\''+\^) 

is selfadjoint. The spectrum of Ma might equal the entire real axis and we note 

Ma{E,H) =i{e~^6H,ij-^dE). 
For more details see [33] or (for the classical case) [32]. 



3 Problem formulation 



We are looking for T-periodic solutions in time of the following generalized Maxwell 
controllability problem 

{dt + iMA){E,H) = {F,G) inS, 

nE = A in T, (3.1) 

{E,H){0) = {Eo,Ho) inn, 

{E,H){T) = iE,Hm inn, 
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where Tt denotes the tangential trace, i.e. Tt 
equation may be written more exphcitly as 



L* in the smooth case. Of course, the first 



dtE-e-'6H 



dtH 



Me 



F 
G 



m ^, 
in S. 



Here / := (0,T) with some time T > is an interval and / = [0,T] denotes its closure. 
Furthermore, we introduce the two product sets S := I x Q and T := / x F. 

As a first order system, the problem at hand represents a natural generalization of 
classical wave equation problems associated to Helmholtz' equation. In [22] Glowinski 
et al. proposed an algorithm to compute time-T periodic solutions u of the prototypical 
scalar linear wave problem 



{d^t -c'A)« = 
Tu = g 
u{0) = u{T) 
dtuiO) = dtuiT) 



m ^, 
on T, 
in Q, 
in Q. 



(3.2) 



They utilized a truncation Qp := |x G Q : \x\ < pjoffi introducing an artificial 
boundary (a sphere Sp of radius p containing ]R^\f2) and a first order absorbing boundary 
condition on it, i.e. setting the translation of Sommerfeld's radiation condition to the time 
dependent formulation {c~^dt + dr)u to zero. Here c is a positive real number and g is 
a given time dependent boundary data. Furthermore, r denotes the usual scalar trace 
operator and r the Euclidean norm on M^. 

They transformed the latter system via the well known substitution 

E := dfU, H := gradu 



into a first order system of 'linear acoustics' 



(9. 



0" 




div" 


1 




grad 



)iE,H) = {0,0) 

tE = dtg 
c-^E + H = 

{E,Hm = {E,H){T) 



in I X fl 



pi 



on T, 
on / X Sp, 
in VLr 



which has a 'Maxwell-type flavor', albeit simpler. Here S,{x) := a;/|a;|. One of the advan- 
tages of this first order system is that it allows for its solution an algorithm using 



L\np) X L\np) 



N 



as control space, i.e. the space of initial data. In former works there was always at least the 
first part of the control space a closed subspace of H^(f2p), which makes the corresponding 
numerics much more difficult due to the need of preconditioning, for instance, in conjugate 
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gradient algorithms. Such preconditioning is not necessary if one uses a purely L^(f2p)- 
control space. 



Utilizing the framework of alternating differential forms, our problem (3.1 ) generalizes 
this approach not only to the classical Maxwell equations in three dimensions but also 
to their generalized and coordinate free version. We should mention that the generalized 
approach also comprises the system of linear acoustics and the 2-dimensional version of 
Maxwell's equations as well as the system of linear elasticity (with another boundary 
condition). 

We emphasize that for g := as well as {F,G) := (0,0), \ := g, e := ^ := 1/c and 



u := E the original problem (3.2) is recovered. 

To start our analysis, we first have to establish a solution theory for the boundary 
value CP 

{dt + iM^){E,H) = {F,G) inS, 

= A in T, (3.3) 

{E,Hm = {Eo,Ho) inn 

with given right hand sides F, G and A as well as initial data [Eq, Hq) belonging to our 
control (Hilbert) space 

H := L^'^'''+'(fi). 



4 Solution theory for the Cauchy problem 



In order to solve (3.3), as a first step we must extend the boundary data from T to S. 



4.1 Traces and extensions 

Recently Week [13] showed how to obtain traces of differential forms on Lipschitz bound- 
aries. Let fib be a bounded Lipschitz domain in with boundary Lb. Then by |l3l 
Theorem 3] there exists a linear and continuous tangential trace operator 

r,,b : D'^in^) D-i/2,,(r^)^ 
where with the notations from |33] 

D-V2,,(r^) 1^ ^ H;iA«(rb) : dr,A G H;iA^+i(rb)}. 

Here dpi, denotes the exterior derivative on Lb. Moreover, by |43i Theorem 4] rt,b is 
surjective and there exists a corresponding linear and continuous tangential extension 
operator (a right inverse of Xt^b) 



't,b 



D-V2,g(r^) _^ D'?((^b). 
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Applying the well known Helmholtz decomposition 

where we introduce the finite dimensional space of Dirichlet forms 

o 

and using rt^bD'^(fib) = {0} we receive a linear and continuous tangential extension oper- 
ator 

f,,b : D-i/^'^(rb) D^in^) n e-'6A'^+\n^). 

Now, we return to our exterior Lipschitz domain Q C M^. Using an usual cut-off-technique 
we obtain the following 

Lemma 4.1 There exists a linear and continuous tangential trace operator 

Tt : D'^{Q) — > D-i/2,g(p) 

and a corresponding linear and continuous tangential extension operator (right inverse) 

ft : D-i/2'«(r) — > D'^iQ) n e-^A^iQ), 

which even maps to forms with fixed compact support and satisfies on D~^/^'^(r) 

Ttft = Id . 

o 

The kernel of equals D'^{Q) and is even well defined on Dl^^{Q). Moreover, ft can 
be chosen, such that supp ftA C ^Ip holds for all X G D^^/^'^(r) and for a fixed p > with 
R^\n cUp. Here, Up C denotes the open Euclidean ball with radius p > centered 
at the origin. 

Remark 4.2 // the boundary is sufficiently smooth, e.g. C"^^^, then even 

holds for all forms E G H™'''^(i7) or E & Moreover, applied to smooth forms from 

C°°''^{Q) we have = l* and, of course, commutates with the exterior derivative. On 
the other hand, if \ e H™-i/2,g(r) we may choose an extension, such that ftA G H™'''^(i7) 
holds and ftA is supported in Vtp. For details see JUJ^ . 

Tt and ft may also be defined on time dependent forms. We get bounded linear 
operators 

rt:F(/, D^(fi)) ^F(/, D-"^'\T)) 
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and 



ft : F(J, D-^/^'^{r)) F(J, D%Q) n e-^A'?(0)) 



with similar properties as mentioned above, where the function space F could be, for 
instance, C^, L^, H^. 

Finally, we also need the corresponding normal trace and extension operators 

r„ := (-1)^^ *r r,*, := (-1)«(^"^) * f,*r 



defined on (g + 1)- respectively (g — l)-forms, where *r denotes Hodge's star operator on 
the (A^ — l)-dimensional submanifold F of Q. 

4.2 Solution theory 



Now we return to the CP (3.3). Let A G C^(/, D ^/^'"^(F)). Then the canonical ansatz 

{E,H) := (E.H) -in\,0) 
leads to a problem with homogeneous boundary condition 



and new data 



{dt + iMj,){E,H) = {F,G) inS, 
nE = in T, 

{E,Hm = {Eo,Ho) inn 



(F, G) := (F, G) + (- dt n\, fi-' d ft A), 
{Eo,Ho) := (Fo,i^o)- (TtA(0),0). 



(4.1) 



Since Ma from (2.1) is linear and selfadjoint, spectral theory suggests a solution {E,H) 



of ( |41| defined for alH G [0, oo) by 

{E,H){t) ■.= exp{-itMA)iEo,Ho)+ I exp {- i{t - s)Ma){F,G){s) ds 



exp(-itMA)( (Fo,^o) + / exp(isMA)(F, G')(s)ds 



Let us analyze this solution thoroughly. For instance, considering forms {Eq, Hq) G EI and 
(F,G) G L2(/,e) we obtain {E,H) in C°(/,e) and thus a solution 



(F,i/)GCO(/,e), 



(4.2) 



if {Eq,Hq) G EI and {F,G) G L^(/, EI). Assuming stronger assumptions on the initial and 
right hand side data, i.e. {Eo,Ho) G D{Ma) and (F, G) G C°(/,e) n (/, F)(Ma)) , we 
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even get a solution {E,H) belonging to C^(/, H) fl C°(J, D(Ma)). Hence, we achieve a 
solution 



if, for instance. 



{F,G) e c°(/,e) n d'{Q) x A^+i(fi)), 

dnXit) e fiA^+^n), tel, (4.3) 
nFit) = dtXit), 
nEo = A(0). 



Then {E, H) is a solution of the CP (3.3) in the strong sense. 
Summing up we obtain: 



Theorem 4.3 Let A G 0{l ,D~^/^^'i{T)) as well as {Eq.Hq) and {F,G) satisfy (Q. 
Then the CP (3.3) is uniquely solved in 

c\i,m) n c°(/, D^iQ) X A''+\n)) 

by 

(E, H){t) := (ftA, 0)(t) + exp(- i ^Ma) (^o - rtA(O), Ho) 

+ / exp (-i(t-s)MA)(F-9,ftA,G + /x"MftA)(s)ds 



/or t G /. We call {E,H) the strong solution of the CP (3.3) with data {F,G, X, Eq, Hq). 

Actually, we are interested in the purely L^-type Hilbert space EI as control space for 
the initial data and even not in /^(Ma) or D''(i7) x A'^^^{Q). Moreover, the constraints 
(4.3) are too complicated and the assumptions on the data much too strong. Thus, we 
have to weaken our solution concept. To approach weak solutions we first have to define 
suitable test forms. 

Definition 4.4 For ($0, ^0) ^ E>{Ma) and t eR the family 

($,^)(t) :=exp(-itMA)(<fo,^o) 
defines a strong solution of the homogeneous Cauchy problem (HCP) 

(9t + iMA)($,^) = (0,0) inRxn, 
rt<l> = mRx T, 

($,^)(0) = (<l>o,^o) ^nn. 

These solutions ($, ^P) are elements of C^(R,M.) fl C°(M, D(Ma)) and we will call them 
test forms with initial values ($0, ^o)- 
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Next, we present the idea of the definition of weak solutions. Thus, let {E,H) be a 



strong solution of (3.3) and be a test form with initial value ($o; ^o) ^ -D(Ma). 

Then we may compute 

((F, G), ($, vl/))^ = {id, + i M^)iE, H), ($, m))^ 

= { dt{E, H), ($, v]/))^ - {M{E, H), ($, ^)>L.,.,+.(^) 
= dt {{E, H), ($, vl/))^ - {(E, H), dti^, ^))^ 
- {dE, - {SH, $)L2..(n)- 

o 

Since $ G D'?(fi) we obtain 

{6H,<l>)L2,^n) = -(if,d$)L2.,+i(n) 

and assuming for these heuristic arguments that E, \1/ and F are sufficiently smooth we 
get by Stokes' theorem 

(dE,$)L2..+i(f^) + (E,5^)L2,5(n) = / d{EA*^) 

Jn 

= j L*{EA = {-ly^ j L*E A *r *r * ^ = (rt^, t^^)l-^ 



2>9(r)- 



Putting all together yields 

((F, G), ($, vl/))^ = 5, ((F, ($, vl/))^ - (A, r,M/)L2.(r) 

-((F,g), (g, + iMA)(«l>,vl/) )^. 



Hence, we only have to remove the time derivative from the forms (F, H) to get our weak 
solutions. (Please compare to Week [12].) 

Definition 4.5 Let (Fq, Hq) E M and (F, G) G L2(/, H) as well as A G (/, D-i/2.9(r)) . 



Then the pair of forms {E,H) is called a weak solution of the CP (3.3) with right hand 







side and initial data {F, G, A, Eq, Hq), if and only if {E, H) belongs to C°(J, H) and 

((F, H), ($, m))^{t) - ((Fo, H,), ($0, v^o)>e 
;(F,G),(<l>,M/)>jj(s) + (A,r,M/)L2.(r)(5)) rfs 

holds for allt E I as well as for all test forms ($, \1') wt/i initial values ($0, ^o) £ F(Ma). 

Remark 4.6 T/ie term (A, rii\I')L2.9(r) needs some detailed interpretation. The normal 
trace of a {q + l)-form from A'^~^^{Q) is only an element of 

^-iA<7(r) := {A G H;i/2'''(r) : 6rX G H-^/^''^-\r)} , 
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where 5r '■= (— l)*-^ ^^^^ *r dr*r denotes the co-derivative on T applied to q-forms. 
Please see again ^431 for details. Hence, at first sight the scalar product 

(A,r,^)L2.,(r)(s) = (A(s),r,^(s))L,,^^^^ (4.4) 

for almost all s makes only sense as an usual dual pairing 

r^^(s)A(s) = (A(s),r„^(s))^i/2..(r)_H-i/2.^(r)- 

Thus, A(s) should be an element of H]/'^''^{T) for almost all s, which is not the case. But, 
since for almost all s the boundary forms A(s) G □^^/^''^(F) and rn\I/(s) G A^-'^/^''^(r) have 



•/tt 

s. We will clarify this in the next lemma. 



more 'regularity' than H^y^ '''(r), the scalar product ( 4.4[ ) still makes sense for almost all 



Lemma 4.7 The {.^''^{T)- scalar product may be extended as a continuous bilinear form 
to D~^/^'''(r) X A^^/^''^(r) (using Stokes' theorem) by the mapping 

b : D-^^^'^iT) X A-^/2,9(r) _^ c 

with 

b{a,f3) = (dfta,fn/3)L2,<,+i(n) + (rta, 5rn/3)L2,,(Q). 
Moreover, for all {E,H) G D'?(fi) x A'^+^{Q) Stokes' theorem 

(dE,ff)L2„+i(n) + {E,6H)^2„^n) = b{T^E,T^H) 
remains valid. Further on we will denote b as usual by {■ , ■ )\_2,q(r)- 

Proof For a G D~^/^'^(r) and (3 G A^-'^/^''^(r) the respective extensions f^a and 
to Q are elements of D'^(f2) and A''+^(r2). Therefore, the definition of b makes sense. 
To show that b is well defined, i.e. does not depend on the extensions, we pick some 
(E, H) G D«(^]) X A5+^(^]) with t^E = a and r^H = (3. Since n{E-na) and T^{H-f^P) 

o o 

vanish we have E — f^a G □'^(fi) and H — fn/3 G A'^^^{Q). Thus, by definition (or an 
approximation argument) 

= ( d{E - fta), ^^>L2,,+i(f^) + {E- fta, 5iJ)L2,,(Q), 
= {dfta,H - fn/3)L2,,+i(Q) + {fta,5{H - rn/3))|_2^^^^ 

hold. Addition shows {dE, H)i2.q+i(^Q^ + {E,6H)\_2,q(^Q^ = b{a,f3), which proves also the 
asserted formula. Finally, the continuity of b follows from the Cauchy-Scharz inequality 
and the continuity of the extensions. □ 



We are ready to prove the main result of this section. 
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Theorem 4.8 There exists at most one weak solution of (3.3). // additionally, for 
instance, A G H^(/, D^^^'^''^{T)^ then there exists always a unique weak solution of (3.3), 
which belongs to C°([0, oo), H) since T is arbitrary. [0, oo) may be replaced by M as well. 

Proof The difference {E, H) of two solutions satisfies 

((E,i7),($,vi/)yt) = o 

for aU t and all test forms ($, Since exp(itMA) is an unitary operator and D(Ma) is 
dense in EI we obtain 

exp{i tM a) {E,H){t) = (0,0) 

and thus {E, H){t) vanishes for all t, which proves uniqueness. To show existence, we use 
the solution {E, H) from Theorem 4.3 suggested by spectral theory, which is still well de- 
fined and still belongs to C°(/, H) by (4.2) even with our weak assumptions. We note that 
we have replaced the stronger constraint A G C^(/, D~^/^''^(r)) by the weaker constraint 
A G Hi(/, D-V2,g(r)) ^ c°(/, D-i/2,g(r)). So, it remains to check if iE,H) satisfies the 
integral equation of Definition 4.5 For this purpose, let ($, \E')(t) = exp(— itMA)($o, ^o), 
t G M, be a test form with ($05 ^0) ^ -D(Ma). We start with the second term in the sum 
of the representation of {E, H) : 



exp(-itMA)(^o - n\iO),Ho),i^,^m^ 
{Eo - nX{0),Ho), ($0, ^0))^ = {{Eo, Ho), ($0, ^o))^ - (^^tA(0), <^o)^..,^^) 
The third term may be handled utilizing Fubini's theorem as follows: 



exp ( - i(t - s)Ma) (F, G) (s) rfs, ($, ^) (t) 

[ { exp(i sMa){F - ds ftA, G + /i-M nX){s), ($0, ^o))^ ds 
Jo 

[ {{F-ds nX, G + fi-'dnX), ($, ^))^{s) ds 
Jo 

iF,G),i<^,^))^{s)ds+ [ {i-d,nX,fx-'dnX),i^,^))^is)ds 

Jo 



We proceed by calculating the last integral. 

- / {dsftX,e^)i2,g^^){s)ds 
Jo 

= -/ ds{ftX,e^)i2,g^Q){s)ds + {ftX,eds^)L2.<i(n){s)ds 
Jo Jo 

= -(ftA,£$)L2,,(n)(t) + (ftA(0),£$o>L2,,(j^) + ^ (ftA,5^)L2,,(n)(s)rfs 
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Hence, we get 



[ {i-dsnX,fi-'dnX),{<^,^))^{s)ds 
Jo 

-(ftA,£$)L2..(n)(t) + (ftA(0),£<l>o>L2,,(n) 
+ / ((ftA,5^)L2,,(n)(s) + (df^A,^) I2,q + 1(^Q){S) ) dS 

Jo ^ V ' 



(A,rj,^)L2,,(r)(s 



by Lemma 4.7 Putting all together completes the proof. 



□ 



4.3 A new notation 

Let us change to a new and shorter notation, which enables us to follow the forthcoming 
arguments and basic ideas more easily. We set := (0, 0) as well as 



n:={E,H), f 

uo := u(0) := (Eo,ifo), ex 
Ut := u(T), gA 



(rtA,0), 



-(9tTtA,yU MftA). 

With this notation our inhomogeneous Cauchy problem (ICP) ( |3.3 ) reads as 



{dt + iM^)n = { 

TtTTU = A 

u(0) = uo 



m ^, 
in T, 
in Q, 



(4.5) 



where for a pair of forms vr denotes the projection onto the first component. Moreover, u 
may be decomposed into u = u' + u'^, where u' and u*^ are the unique weak solutions of 
the CPs 



(dt + iMAW = 0, 

TtW = 0, 

u'(0) = Uo, 



(at + iMA)u^ = f 

TtTTU'^ = A 

u'=(0) = 



m ^, 
in T, 
in Q. 



(4.6) 



u' depends linearly and continuously on the initial data uq and u'^ is independent of the 



initial data uq, i.e. constant with respect to Uq. The unique weak solutions of (4.5) and 
( [46| exist by Theorem IZsl in C°(7,H) for all T and all 



Uo G e, f G e), A G (/, D-^/^'%T)) 



(4.7) 
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and are given by the following formulas: 



u(t) = e,(t) + e-'™- (uo - e,(0)) + + g,)(,) ds 

Jo 



u'(t) =e-'™'^ Uo 

Jo 



(4. 



■ i{t-s)MA 



(f + gA)(s) 



5 Least-squares formulation of the controllability prob- 
lem 



From now on, let the right hand side data f and A satisfy (4.7) as well as the time 
T > be given and fixed. 



In order to solve the controllability problem (3.1), which reads now 



'Find uq G H, such that u satisfies (4.5) and = uq. 
we investigate the equation 

— Uo = 



more thoroughly. With the help of (4.8|) we obtain 

= 

-^™'^-l)uo + u^ 



ut = u(T) = u'(T) + u%T) = e-'™^ uq + u^. 



(5.1) 
(5.2) 



Ut - Uo = 

Consequently, with the continuous linear operator in EI 

Ct ■.= C{t) :=e-'*^^^-l, 
which satisfies ||Ct|| < 2 for all t and will be called 'control operator', we get 

Ut - Uq = CtUo + Uy. 

Hence, we have to solve the linear equation 

CtUq + Ut = 



in the Hilbert space H, which we want to try approximately by an CGA. Since, of course, 
Ct is neither symmetric nor selfadjoint the usual CGA suggests to consider the corre- 
sponding normal equation 



(5.3) 



CtCtUq + CtUt — 



(5.4) 
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where Q = e'™^ -1 denotes the adjoint operator of Ct- We note Q* = Ct and that C^Ct 
is selfadjoint. Consequently, we are forced to consider and to minimize the quadratic 
functional J-" with 

J'(uo) := ^(QCtUo,Uo)h + 3fJ(CyU^,Uo)H 

which, of course, is minimized, if and only if the quadratic functional T := T ^ ||u^||^/2, 

7 : W [0,00) 

I > \\n „ I ,,c ||2 _ 111., 112 5 W-^J 

is minimized. This leads to the following least-squares formulation: 
'Find initial data uq G H, such that 

VvoGH ^(uo)<^(vo).' (5.6) 



Here, u respectively v is the unique weak solution of the ICP (4.5) with initial data uq 
respectively vq. 

The implementation of the CGA in EI is greatly facilitated by the knowledge of the 
derivative T' . Since T is differentiable as a quadratic functional we get from 

J-(uo + vo) = -F(uo) + 3f?(CTVo,CTUo + u^)e + ]^CT^^f^. (5.7) 
where uo,vo G H, immediately 

^'(uo)vo = 3?(vo, C* (CtUo + v5^))^ = 3ft (vq, C^C^uo + C*u^)h (5.8) 
and, of course, the normal equation is recovered. In this sense, we may identify 

J^'(uo) with C^CtUo + Qu^ G M. 
Furthermore, we receive the representations 

Vt := V{t) := C;Ct = (e'*^'^ -l)(e-'™^ -1) = 2(l - cos(tMA)), 

:= u(t) := C*u^(t) (5.9) 

= (ei*MA -l)ex{t) + (e-'™^ -l)eA(O) + [\l - e"^™^) e'^^^f + gx){s) ds, 

Jo 

where we will call the continuous linear operator Vf in EI the 'derivative operator'. We 
have W'DtW < 4 for all t. Finally we obtain 

J"(uo)vo = 3ft(vo, VtUq + ut)r. (5.10) 



By (5.7) we get also 

^(uo) < ^(uo + Vo) - J"(uo)vo 

for all Vq G EI and thus 
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Remark 5.1 For uq G EI the following assertions are equivalent: 

(i) Uq is a solution of the least squares problem (5.6). 

(ii) J-'(uo) = 

(iii) T^tUq + ut = [normal equation (5.4)) 



Using (5.3), let us interpret the derivative vector 
VtUo 



Ut = C^Ug G H, Ug := - Uo G EI 
more thoroughly. Clearly, the forms u*'"*" and u*'~ defined by 

u*'+(t) := e'*^^ u*, u*-(t) := e'^^-*)^^ u* 
are the unique weak solutions of the homogeneous adjoint Cauchy problems (HACPsi) 





ul, u*'-(T) 



(9iTiMA)u*'± 



u*'+(0) 



u 







m ^, 
in T, 
in Q 



(5.11) 



and we have u*'+(T) = u*'"(0) = e 



u5. I.e. 



u 



T 







u 







u 



0- 



Here, the signs ± indicate that the wave u*'+ evolves forward in time, whereas the wave 
u*'~ evolves backward in time. Of course, this implies a change of sign in the dt-teTm. 



We note that we define the weak solutions of the HACPs± analogously to Definition 4.5 



Finally, we obtain two more nice representations of our derivative vector utilizing the 
solutions of the HACPsi ( |5TI| 

"Dj-Uq + = — Uq = Uq ^ — Uq. (5.12) 

As already pointed out, the derivative vector depends on the initial condition uq both 
directly and indirectly through the solution u of the ICP (4.5) and one of the solutions 
u*'^ of the HACPsi (5.11). Moreover, we saw in (4.8) that u = u' + u'^ splits up into a 
linear and continuous and a constant part (with respect to uq). Of course, the same holds 
true for the solutions of the adjoint equations. Let us pick, for instance, the forward in 
time solution u* := u*'"*". Then u* depends linearly and continuously on the initial data 
Uq and may be decomposed into u* = u*'' + u*'^, where u*'' and u*'^ are the unique weak 
solutions of the HCPs 



(9t-iMA)u*''/^ 

*,//c 



Ttvru 



(0) 






u 



*,l/c 



m ^, 
in T, 
in Q 
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with Uq' := u(p — uo = CtUq and Uq'^ := as well as Uq = Uq'^ + Uq'^. Again, u*'' depends 
linearly and continuously on uq, whereas u*''^ does not depend on uq. Of course, we have 

Putting all together, we see 

Z^TUq = — Uj. + Uq, Ur = — Urp. 

6 Conjugate gradient algorithm for the least-squares 
problem 



Although it has become customary to use CGAs in Hilbert spaces, see e.g. [HI [TTl 
\TE\ HSl [201 Ell [22] as a selection, we briefly want to repeat the algorithm here. 

In order to solve approximately our least squares problem (LSP) ( |5.6 ), i.e. the linear 
equation 

CtUq + = 



or by Remark 5.1 equivalently our normal equation (5.4) 

VtVLq + ut = C^CtUo + C^u^ = 

we will use the following variant of the usual CGA: Given an approximation Uq~^ and 
last search direction d"~^ we compute the new search direction and approximation by 



„n-l 



n— 1 jn— 1 



U 



•" 



u 



n-1 



with coefficients 



a 



/3" : = 



l|CTd"||^ 

WCTd-Wl 



„n-l||2 



l|CTd"||^' 



„n\\2 



.nil 2 



||/^„f4n II 2 ||„n— 1||2 ' 

a \\Lt(1 lie F \\m 
where the residual is given by 

r := C^-CtUq + C^u^^ = r + a C^UTd . 

We note that the initializing procedure in the CGA 



(6.1) 



I.e. 



r*n* 



u 



T 



U 



0' 



Un 
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where we picked the forward in time solution u* := u*'"*", needs the solution u at time T 



of the ICP (4.5) with initial data Uq as well as the forward in time solution u* at time T 
of the HACP+ (5.11) with initial data Uq. Analogously the procedure within the loop of 
the CGA 



d = C^Crd", 



(6.2) 



I.e. 



d = C^ul 



u 



T 



0' 



Cyd" = ut- d", 



where we again used the forward in time solution u* := u*'+, needs the solution u := u' 



at time T of the HOP (4.6) with initial data d" as well as the forward in time solution u* 



at time T of the HACP+ (5.11) with initial data u* 



We recall that the procedure (6.1) respectively (6.2) may be identified with the calcu- 



lation of the derivative or 'gradient' 



of the least squares functional 



^:=^u^ : H [0,oo) 

uo I — y |||CTUo + u^||e 



respectively 



^0 



[0,oo) 

2ll^TUo||e 



We will present the CGA for the approximate solution of the LSP as our Algorithm [T] 
In the beginning of the algorithm, before entering the iteration loop, we choose an initial 
control vector Ug G EI and compute the first residual vector r'', i.e. the 'gradient' of the 
functional J-'u^ at the point Uq, which gives the first minimizing direction d^ = r*'. The 



computation of this residual requires the solutions of the ICP (4.5) with initial control 



vector Uq and of the HACP+ (5.11). Then, on each CGA iteration we calculate the 



solutions of the HCP (4.6) with initial vector d" and of the HACP+ (5.11). This gives 



the 'gradient' of the functional J-'o at the point d", which is needed to update the new 
residual vector r" and the new control vector Ug. Finally we set the new minimizing 
direction d"^^. 
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Algorithm 1 CGA in H for LSP §^ 



initialization 

set n = 

set initial control vector Uq G EI 



solve ICP (4.5) with initial vector Uq and get u 

solve HACP+ (5.11) with initial vector Uq = ut — Ug and get u* 

compute residual vector (gradient J-'^c^(uq)) r" = u^ — Uq 

compute norm = ||r"||g 

if p" small then 

goto exit 
end if 

set first minimizing direction d*^"*"^ = r" 

loop {for n > 1 assuming Ug~^ and r"^^ 7^ 0, p"^^ and d" 7^ are known} 
solve HCP (4.6) with initial vector d" and get u 
solve HACP+ ( 5.11| ) with initial vector Ug = ut - d" and get u* 
compute gradient (J-'o(d"')) d = u^ — u* 



-P" 



Un 



„n\\2 



compute parameter a = 
update control vector Ug 
update residual vector r" 
compute norm p" = ||r 
if p" small or n large then 

goto exit 
end if 

compute parameter p = l/p""^ 
compute parameter p = p"p 
update minimizing direction d"+^ 
set n = n + 1 

end loop 

exit 

take Ug as solution 



V(d,d")i 



ad" 
ad 



r" + pd" 



We note that we may use the backward in time system HACP- (5.11) instead of 
HACP+ as well. Then, in this variant by (5.12) we have to replace the computation of 
the residual or gradient vector u^ — Ug = u^ — Ug by Uq'~ — Ug. 



7 Translation to classical problems 



We briefly mention, which classical problems of vector analysis are covered by our 
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general CP (3.3). (3.3) reads: 

dtE - e-^6H = F in S 

dtH - fi-^dE = G inS (7.1) 

nE = A in T 

The initial condition always stays the same. 

In the exterior derivative d and the co-derivative S turn to the classical differential 
operators from vector analysis 

grad = V, curl = V x , div = V • 

and the well known standard Sobolev spaces 

l^{n), H(grad,n) = Hi(fi), H(curl,fi), H (div, fi) 

appear. Moreover, the tangential trace becomes the usual scalar, tangential or normal 
trace, respectively. As long as the operators x or curl are not involved, the classical 
calculus extends to M^, e N. 

We obtain the following problems in : 

• q = N (trivial case) : dtE = F 

• q = (linear acoustics, Dirichlet case): 

dtE - e'^ divH = F in S 

dtH — grad E = G in S 

E\-p = A on T 

• g = — 1 (linear acoustics, Neumann case): 

dt E — grad H = F in S 

dt H - fi-^ div E = a inS 
u ■ E\j. = A on T 

• q = 1 and A^ = 3 (Maxwell's equations): 

dtE + e-^cm\H = F in S 

dtH — /i"^ curl E = G in S 

u X E\-r = A on T 



We note that the equations of linear elasticity are also covered by our approach, if 
we change the tangential boundary condition into the more simple one of componentwise 
scalar Dirichlet boundary conditions. See, for example. Week and Witsch jSj. 
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8 Conclusion and outlook 



The considered approach of combining the exact controllabihty method with DEC- 
based discretization seems to be a promising way to compute time-harmonic scattered 
waves. The problem setup is general enough to treat the most important cases of lin- 
ear wave propagation, i.e. electro-magnetic, acoustic and elastic waves in three space 
dimensions. 

There are certain key benefits of the method. The DEC approach leads to a discrete 
scheme that has good conservation properties [16]. Also, the resulting time integration 
scheme can be implemented in explicit manner. As the mass matrices are also diagonal, 
the time integration will be computationally very efficient and all the related computations 
are easily parallelized by standard domain decomposition techniques with coarse grained 
boundary swapping approach. The parallel implementation of the outer CG-iterations 
is also very straightforward because we minimize the error in periodicity in the squared 
energy norm of the system. In the family of problems we consider, the energy norm is 
a weighted L^-norm, which means that the discrete quadratic functional we minimize is 
spanned by a diagonal mass matrix. Hence, in practice, no preconditioning is required. 
This is supported by our initial experiments (appendix), at least when the mesh is refined. 
The convergence of CG iterations did not depend on the mesh step size. 

The use of DEC causes some challenges to overcome. These are related to the defini- 
tion of the dual mesh and the resulting discrete Hodge operator. In our initial experiments 
(appendix), we used well-centered meshes and the circumcenter based dual mesh defini- 
tion, which naturally leads to a diagonal discrete Hodge operator. The tiling of a general 
domain with well-centered simplices is an open problem, which has been accomplished for 
some simple shaped domains [H]. 

To make the method applicable to general cases arising from practical problems, we 
must allow for simplicial meshes, which are not well-centered. Therefore, e.g., the barycen- 
tric dual meshes need to be considered. It should be noted that simplicial meshes are not 
obligatory in our approach. As was shown in \iQl, the classical and very widely used 
Yee-scheme for Maxwell's equations |16] is just a special case of general DEC-approach 
and the same control scheme can be implemented directly for Yee's scheme as well. 

A Appendix: Preliminary numerical results 

We have implemented the DEC for 3-dimensional geometries to solve electro-magnetic 
problems. Following [16l section 4.3] the implementation allows us to have an unstruc- 
tured mesh in space and asynchronous time steps. Our implementation is based on the 
circumcentric dual mesh, which is the most simple way to build a DEC solver, but it 
also requires Delaunay's property of the mesh. Since we are interested in time-periodic 
solutions, we implemented the CGA using the theory of section [6] as well. 
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A: 

— 1^ 1^^^ — 1 — 


— 1 — 1 — 1 — 1 — i-^ 






\ 




Electric 


field E 



yi 




Magnetic field B 



Figure 1: We see simulated fields at t = T. The pictures are cross sections on a xy-plane, 
where the x-, y- and z-components of the field vectors are presented in blue, green and 
red, respectively. The zero field would be 50% gray. 



We discuss some preliminary results of our simulations. Let us consider a scattering 
problem, where electro-magnetic plane waves are reflected by a sphere and scattered to 
infinity. We are interested in the accuracy of the simulation and in the convergence of the 
CGA. 




number of iterations 



Figure 2: convergence of the residual in CGA 
For our very simple model radiation problem 



curli/-ia;^ = 
curlE" + iuH = 

u X E\r. = A 



in Q, 
in Q, 

on r, 
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Figure 3: relative error of the simulated electric field at time t = T integrated over the 
mesh volume 

where Q := |a; G : |a;| > l} is the exterior of the closed unit ball, i.e., with a 
spherical hole (ball) of radius 1 in the middle, T = S"^ and where we picked the frequency 
u = 27r/3, the exact solution is known explicitly and can be found, for instance, in the 
book [13, Theorem 6.25]. We took only one non-zero component = 1. This simplifies 
the solution to 

H = - curl E, E = curl E, E = h{ujr)y{^) Id 

CO 

with Hankel's spherical function of first kind h := h^i^ and spherical harmonic y := y^. 
More explicitly, we have 

= /i(a;r)r(0 X e, Y ■.= Vy, 

where V denotes the spherical gradient on F. Let us note that 

lyx E = h{uor)i x ^ x Y{i) = h{ur)Y{^) 

since u = and Y{^) is tangential at F. We generate a wave on the boundary F by 
the Dirichlet boundary condition A := h{u)Y, i.e., setting r := 1. Picking an artificial 
outer boundary F, the sphere of radius 5 centered at the origin, we impose the classical 
Silver-Miiller first order absorbing boundary condition 

E + H = out. 

We have simulated the test problem with six different meshes of varying element sizes. 
The initialized edge lengths varied from about 1/5 to 1/2. In Figure |2] we see how the 
residual converges in the CGA. After 140 loops of the CGA we got about 10"^ times 
smaller residuals. The convergence seems to be independent of the mesh element size. 
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In Figure [3] we plotted the differences between the simulation corresponding to different 
mesh element sizes and the exact solution. The error of the simulated fields is decreasing 
when the mesh is refined. The decrease of the error even might be of second order with 
respect to the average edge length. 
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